*------------------------------------------------------------------------------*
* This program analyzes simulated data frz at 58 with average HRS replacement
*------------------------------------------------------------------------------*

global indata ".../data"
global sigma = 1.1342 //estimated EIS
global r = .029 //real interest rate
global zage = 58 //freeze age
global dage = 60 //age for discounted PVs

*Simulated data
import delimited using "$indata/simulated_profiles_z58_avg.csv", clear

*Labels 
rename v1 persid
rename v2 g
rename v3 emp
rename v4 A
rename v5 dca
rename v6 y
rename v7 db_ben
rename v8 ssben
rename v9 mortp
rename v10 trt_flag

*Current age
gen sortord = _n
sort persid sortord trt_flag
bysort persid trt_flag: gen n = _n
gen t = n-6
gen age = t+58
drop n t

*Psuedo persid
egen persid2 = group(persid trt_flag)
xtset persid2 age

*pension accruals
gen dca_1 = l.dca
gen dcaccrual = dca - dca_1
gen dcaccrual_f = f.dcaccrual

*deferral fraction
gen defrate = dcaccrual_f/y
gen mw = (defrate-.011)/1.39 if defrate<.10135
replace mw = (defrate - .03635) if defrate>=.10135
replace mw = 0 if emp==0
replace mw = max(0,mw)

*Deferral amount
gen ern_after_dc = (1-mw)*y

*Merge with ann factors
merge m:1 age using "$indata/sim_annfactor", keepus(annfactor)
keep if _merge!=2
drop _merge

*Earnings defined only when employed
replace ern_after_dc = 0 if emp==0

*Liquidated DC wealth 
gen dcdraw = -1*dcaccrual_f
replace dcdraw = 0 if emp==1
replace dcdraw = 0 if dcdraw==.
*Interpolation error sometimes causes dcdraw to be <0 in a small # of cases
replace dcdraw = max(0,dcdraw)

*Retirement age
gen temp = age*emp
bysort persid2: egen retage = max(temp) //for each person
drop temp
*Ret age is the first age out of LF
replace retage=retage+1

*ss ben at retirement 
gen ssbenR = ssben if age==retage //for workers who retire before 62
gen ssben62 = ssben if age==62
bysort persid2: egen temp = max(ssbenR)
replace ssbenR = temp if ssbenR==. & emp==0 & retage>=62
drop temp
bysort persid2: egen temp = max(ssben62)
replace ssbenR = temp if ssbenR==. & emp==0 & retage<62
drop temp
replace ssbenR = 0 if emp==1
replace ssbenR = 0 if age<62
drop ssben62

*db ben at retirement
gen db_benR = db_ben if age==retage 
bysort persid2: egen temp = max(db_benR)
replace db_benR = temp if db_benR==. & emp==0 
drop temp
replace db_benR = 0 if emp==1

*File to pass to TAXSIM to compute FIITAX liability

//Working individuals 
preserve

	keep if emp==1 
	keep sortord age ern_after_dc 
	
	*file structure for taxsim

	gen taxsimid = sortord   
	gen year = 2010     //2010
	gen state  = 0     
	gen mstat = 1       //1 = single 
	gen page = .        //age
	gen sage = 0       
	gen depx = 0
	gen dep13 = 0
	gen dep17 = 0
	gen dep18 = 0
	gen pwages = .      //annual labor earnings
	gen swages = 0
	gen dividends = 0
	gen intrec = 0
	gen stcg = 0 
	gen ltcg = 0
	gen otherprop = 0
	gen nonprop = 0     
	gen pensions = 0    //pension distributions
	gen gssi = 0        //SS benefits
	gen ui = 0
	gen transfers = 0
	gen rentpaid = 0
	gen proptax = 0
	gen otheritem = 0
	gen childcare = 0
	gen mortgage = 0
	
	replace pwages = ern_after_dc
	replace page = age
	drop age ern_after_dc
	
	*Run the tax calculator
	cd "$indata"
	taxsim35, replace

	*Keep relevant vars
	keep fiitax fica sortord
	save "$indata/z58_avg_Wsiminput_fiitax", replace
restore

//Retired individuals 
preserve

	keep if emp==0 
	keep sortord age dcdraw db_benR ssbenR
	
	*file structure for taxsim

	gen taxsimid = sortord   
	gen year = 2010     //2010
	gen state  = 0     
	gen mstat = 1       //1 = single 
	gen page = .        //age
	gen sage = 0       
	gen depx = 0
	gen dep13 = 0
	gen dep17 = 0
	gen dep18 = 0
	gen pwages = 0      //annual labor earnings
	gen swages = 0
	gen dividends = 0
	gen intrec = 0
	gen stcg = 0 
	gen ltcg = 0
	gen otherprop = 0
	gen nonprop = 0
	gen pensions = .    //pension distributions
	gen gssi = .        //SS benefits
	gen ui = 0
	gen transfers = 0
	gen rentpaid = 0
	gen proptax = 0
	gen otheritem = 0
	gen childcare = 0
	gen mortgage = 0

	replace pensions = dcdraw + db_ben
	replace gssi = ssben
	replace page = age
	drop ssben dcdraw db_ben

	*Run the tax calculator
	cd "$indata"
	taxsim35, replace

	*Keep relevant vars
	keep fiitax fica sortord	
	save "$indata/z58_avg_Rsiminput_fiitax", replace
restore


*Merge tax data back 
merge 1:1 sortord using "$indata/z58_avg_Wsiminput_fiitax", keepus(fiitax fica)
drop _merge
rename fiitax fiitax1 
rename fica fica1
merge 1:1 sortord using "$indata/z58_avg_Rsiminput_fiitax", keepus(fiitax fica)
drop _merge
rename fiitax fiitax2
rename fica fica2 

forvalues i = 1/2{
	replace fiitax`i' = 0 if fiitax`i'==.
	replace fica`i' = 0 if fica`i'==.
}

egen fiitax = rowmax(fiitax1 - fiitax2)
egen fica = rowmax(fica1 - fica2)
drop fiitax1 fiitax2
drop fica1 fica2 

*Comparative statics: 
* identify those who retire earlier than no-freeze 
* identify those who retire at the same age as no-freeze
* identify those who retire later than no-freeze
* then compute \delta retage by persid 

sort persid age trt_flag

gen d_retage = retage - retage[_n-1] if persid==persid[_n-1] & age==age[_n-1]
gen d_retage_flag = 0
replace d_retage_flag = 1 if d_retage>0
replace d_retage_flag = -1 if d_retage<0
replace d_retage_flag = . if trt_flag==0
tab d_retage_flag if trt_flag==1 & age==$zage //one obs per person 

*retirement horizon: (retage - d_retage) - freeze age
gen ret_horizon = retage - d_retage - $zage
tabstat ret_horizon if trt_flag==1 & age==$zage, by(d_retage_flag)

*------------------
*Compute EV
*------------------

*Consumption
gen worker_tax_liab = 0
replace worker_tax_liab = fiitax + fica/2
xtset persid2 age
gen A_1 = l.A
gen deltaA = A_1-A
gen cons = 0
replace cons = (ern_after_dc - worker_tax_liab) + deltaA if emp==1
replace cons = (dcdraw - fiitax) + db_benR + ssbenR + deltaA if emp==0
replace cons = 0 if cons<0 //small number of cases where interpolation makes C <0 
gen disc_cons = cons*(1-mortp)*(1+$r)^($dage-age) //
bysort persid2: egen PVcons = sum(disc_cons) if age>=$zage 
*difference in consumption PV
sort persid age trt_flag 
gen d_PVcons = PVcons - PVcons[_n-1] if persid==persid[_n-1] & age==age[_n-1] & age>=$zage

*Utility
gen utility = 0
replace utility = cons^(1-($sigma)^-1)/(1-($sigma)^-1) - g if emp==1
replace utility = cons^(1-($sigma)^-1)/(1-($sigma)^-1) if emp==0

*EV is income equivalent of d_utility
sort persid age trt_flag
gen income_equiv = (utility*(1-($sigma)^-1))^(1/(1-($sigma)^-1))  //income equivalent of utility in each state 
replace income_equiv = 0 if income_equiv ==.
gen ev = income_equiv - income_equiv[_n-1] if persid==persid[_n-1] & age==age[_n-1] & age>=$zage

*Present value of EV diffs
gen disc_ev = ev*(1-mortp)*(1+$r)^($dage-age)
bysort persid2: egen PV_ev = sum(disc_ev) if age>=$zage 

*------------------
* DC + non-pension wealth at freeze age
*------------------

gen tot_wealth = A+dca //liquid wealth 

*------------------
*studentize g in the year of the freeze
*------------------

gen stud_g = .
summ g if age==$zage & trt_flag==1, d 
replace stud_g = (g - r(mean))/r(sd) if age==$zage & trt_flag==1

*------------------
*Table to output
*------------------

tabstat stud_g tot_wealth PV_ev retage d_retage if trt_flag==1 & age==$zage, by(d_retage_flag)


matrix D = J(4,3,.)

*shares by type
count if trt_flag==1 & age==$zage
local N = r(N) 

*early retirement
count if d_retage_flag==-1 & trt_flag==1 & age==$zage
matrix D[1,1] = r(N)/`N'

local vars "d_retage PV_ev"
local i = 2
foreach v of local vars{

	if `i' == 2{
		summ `v' if trt_flag==1 & age==$zage & d_retage_flag==-1
		matrix D[1,`i'] = r(mean)
	}
	
	else if `i' == 3{
		summ `v' if trt_flag==1 & age==$zage & d_retage_flag==-1, d
		matrix D[1,`i'] = r(p50)
	}

	local i = `i'+1	
} 

*no change in retirement age
count if d_retage_flag==0 & trt_flag==1 & age==$zage
matrix D[2,1] = r(N)/`N'

local vars "d_retage PV_ev"
local i = 2
foreach v of local vars{

	if `i' == 2{
		summ `v' if trt_flag==1 & age==$zage & d_retage_flag==0
		matrix D[2,`i'] = r(mean)
	}
	
	else if `i' == 3{
		summ `v' if trt_flag==1 & age==$zage & d_retage_flag==0, d
		matrix D[2,`i'] = r(p50)
	}

	local i = `i'+1	
} 


*delayed retirement age
count if d_retage_flag==1 & trt_flag==1 & age==$zage
matrix D[3,1] = r(N)/`N'

local vars "d_retage PV_ev"
local i = 2
foreach v of local vars{

	if `i' == 2{
		summ `v' if trt_flag==1 & age==$zage & d_retage_flag==1
		matrix D[3,`i'] = r(mean)
	}
	
	else if `i' == 3{
		summ `v' if trt_flag==1 & age==$zage & d_retage_flag==1, d
		matrix D[3,`i'] = r(p50)
	}

	local i = `i'+1	
} 

*overall average
matrix D[4,1] = 1

local vars "d_retage PV_ev"
local i = 2
foreach v of local vars{

	if `i' == 2{
		summ `v' if trt_flag==1 & age==$zage 
		matrix D[4,`i'] = r(mean)
	}
	
	else if `i' == 3{
		summ `v' if trt_flag==1 & age==$zage, d
		matrix D[4,`i'] = r(p50)
	}

	local i = `i'+1	
} 

preserve
	clear
	svmat D
	gen var = ""
	replace var = "Early retirement" if _n==1
	replace var = "No change in retirement" if _n==2
	replace var = "Delayed retirement " if _n==3
	replace var = "Average" if _n==4
	order var D1
	rename D1 share
	rename D2 delta_R
	rename D3 ev
	export delimited ".../model-output/z58_avg_cf.csv", replace
restore
